home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 441 / fplib20 / trig.c < prev   
C/C++ Source or Header  |  1990-11-23  |  2KB  |  80 lines
  1. /* SIN, COS, and TAN functions for Sozobon C.                */
  2. /* Copyright © David Brooks, 1989 All Rights Reserved            */
  3.  
  4. #include <fplib.h>
  5. #include <errno.h>
  6.  
  7. /* SIN:
  8. /*                                    */
  9. /* The method is from Hart et al: "Computer Approximations". It reduces    */
  10. /* the range to 0..pi/2, scaled to 0..1, and applies:            */
  11. /* 1.57079629t - 0.64596336t^3 + 0.0796884805t^5 - 0.0046722279t^7    */
  12. /*        + 0.00015082056t^9 (precision of 8.5 digits) and    */
  13. /* finally fixes the sign.                        */
  14. /*                                    */
  15. /* MODIFICATIONS:                            */
  16. /* DWB    03/15/89    Don't try to negate x if it's zero.        */
  17. /* DWB    04/30/89    Fix cos(0.0)                    */
  18.  
  19. float sin(a) fstruct a;
  20. {    register float        t, ipart, tsq;
  21.     fstruct            x;
  22.     register unsigned long    tint;
  23.     register char        sign;
  24.  
  25.     sign = a.sc[3] & 0x80;
  26.     a.sc[3] ^= sign;        /* sin(negative x) = -sin(-x) */
  27.  
  28. /* If the argument is huge, the result is rather arbitrary -- and,    */
  29. /* besides, the following code will break.                */
  30.  
  31.     if (a.sc[3] >= (BIAS + MANT_BITS))
  32.     {    errno = ERANGE;
  33.         return 0.0;
  34.     }
  35.  
  36.     t = a.f * M_2_PI;        /* Express in quarter-turns */
  37.  
  38. /* Get integer part and fraction.  Optimise the common case, in the
  39.    first quadrant */
  40.  
  41.     if (t <= 1.0)
  42.         tint = 0L;
  43.     else
  44.     {    tint = t;            /* tint = whole quarterturns */
  45.         ipart = (float)tint;
  46.         t = (tint & 1)?ipart - t + 1.0:t - ipart;   /* Get fraction */
  47.     }
  48.     tsq = t * t;
  49.     x.f = t * (1.57079629 + tsq * (-0.64596336 + tsq * (0.0796884805
  50.         + tsq * (-0.46722279e-2 + tsq * 0.15082056e-3))));
  51.  
  52.     if (x.sc[3] != 0)            /* Don't negate zero!    */
  53.     {    if ((int)tint & 2)        /* Quadrants 3 and 4    */
  54.             x.sc[3] |= 0x80;    /* Set negative        */
  55.         if (sign != 0)
  56.             x.sc[3] ^= 0x80;    /* Negate        */
  57.     }
  58.  
  59.     return x.f;
  60. }
  61.  
  62. /* COS:                                    */
  63. /*                                    */
  64. /* Although approximations are available, this will do fine.  I'm not    */
  65. /* proud of the hack to fix the only rational argument with a rational    */
  66. /* result, but then pride is the first deadly sin...            */
  67.  
  68. float cos(a) float a;
  69. {
  70.     return a==0.0?1.0:sin(M_PI_2 - a);
  71. }
  72.  
  73. /* TAN:                                    */
  74. /* (v-e-r-y s-l-o-w-l-y)                        */
  75.  
  76. float tan(a) float a;
  77. {
  78.     return sin(a) / sin(M_PI_2 - a);
  79. }
  80.